# fits the two-stage model which is also used to initialise the ego-ERGM EM algorithm.


init.egoergm<-function(ego.terms, G, p.ego.terms=NULL) #Specify function in terms of ego.terms and G
  {
  Nterms<-length(ego.terms) #specifying a "nterms" object as the length of the number of ego.terms
  ergmformula <- paste("~", paste(ego.terms,collapse="+"),sep="") 
      #Specify object ergmformula - paste command has three arguments - the stuff you want to paste together
      #, sep is how to separate them, and collapse is if you want smush it together.  This specifies ergmformula
      #as ~ pasted together with ego terms, separated by " " 
    
  form<-ergm.update.formula(as.formula(paste("x[[i]]",ergmformula)),x[[i]] ~ .)
      #This specifies object "Form" as an object that is a new formula that is a function of an indexed xi, according to 
      #the ergm formula that is specied above
  
  ######### fit all ego-networks #############
  if (is.null(p.ego.terms)) #if p.ego terms is null or empty, then do the following...
    {
    theta<-matrix(0,N,Nterms) #Theta is set equal to a matrix filled with zeros, with N rows and Nterms columns
    for (i in 1:N) #For loop indexed from 1 to N - for every indexed i in the theta matrix, do ergmMPLE in the form object,
      #with the output of the fit, which is a coefficient for the model terms for each N.  
      theta[i,]<-ergmMPLE(form,output="fit")$coef
    # next couple of lines very ad-hoc but not an issue post EM convergence.
    theta[is.na(theta)]<-0
    theta[theta==-Inf]<- -1e6
    theta[theta==Inf]<- 1e6
    ############# initialisation ##############
    # k-means Clustering start #if statement - if G>1, then the following is true and initial.clusters is set equal to
    #the kmeans of which takes the theta matrixm, for the number of centers, it is a function of theta, 2, and some other stuff?
    #I don't really understand this portion too much - go back to the article and see where this comes in? 
    if (G>1)
      {
      initial.clusters<-kmeans(theta, G, centers=apply(theta, 2, tapply,cutree(hclust(dist(theta)),G), mean),nstart=5)
      group.theta<-initial.clusters$centers
      group.theta<-matrix(group.theta,G,Nterms)
      }
    if (G==1)
      group.theta<-matrix(apply(theta,2,mean),nrow=1)
    lambda<-matrix(NaN,N,G)
    for (i in 1:N)
      for (g in 1:G)
        lambda[i,g]<-1/(sqrt(t(group.theta[g,]-theta[i,])%*%(group.theta[g,]-theta[i,]))+tol) # nugget to avoid zero
    lambda<-lambda/apply(lambda,1,sum) # normalise lambda
    print("Finished kmeans initialisation")
    }
  return(list(theta=theta, group.theta=group.theta, lambda=lambda, G=G))
  }
